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SUMMARY 

This paper treats the RMS-errors associated with the position and ve- 
locity of a satellite or spacecraft when tracked by all types of present day 
tracking systems. These errors are based on uncertainties in measurements 
made with the systems as well as those associated with their location. 

The present paper (Part I) is principally a theoretical treatment which 
establishes the mathematical models necessary to solve for the errors in 
satellite position and velocity. It is presumed throughout the paper that these 
errors are to be determined for discrete points in the satellite's orbit. This 
approach enables one to calculate the error propagation during short time 
intervals (order of seconds). This is of particular interest for instance for 
evaluation of a guidance system during a short burning phase. A least square 
solution of non- simultaneous observations would diverge (matrices involved 
become ill-conditioned). This condition imposes a constraint on the method 
of solution which is, that either one tracking system can measure both posi- 
tion and/or velocity, or that several tracking systems observe the satellite 
simultaneously to produce the equivalent effect. Both these alternatives are 
considered in Part I . 

A rigorous derivation of the ' Method of Least Squares'' is also presented 
for completeness, since it is used to a large extent. 

Part II (presently in preparation) will show in detail the application of 
the equations derived in Part I assuming simultaneity as well as non- 
simultaneity. In the latter case position and velocity need not be fully deter- 
mined by the tracking system or system complexes. A number (i>6) of 
range, range- rate, or angular measurements are adequate. Thus the con- 
straint mentioned above and applicable for Part I does not exist, hence mak- 
ing the method more general. Numerical examples and results based upon 
known errors associated with the systems and systems locations will be pre- 
sented for simultaneous and non-simultaneous derivations. 
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TRACKING SYSTEMS, THEIR MATHEMATICAL 
MODELS AND THEIR ERRORS 
PART I - THEORY 


by 

F. O. Vonbun and W. D. Kahn 

Goddard Space Flight Center 


INTRODUCTION 

In recent years, requirements of satellite and spacecraft tracking accuracies have 
increased considerably; consequently more and more tracking systems and stations have 
been employed. From this situation, the question evolved: How accurately can the posi- 
tion and velocity of a space vehicle to be determined when it is tracked by a combination 

of tracking systems? 

This paper attempts to present a detailed study of the propagation of errors in data 
obtained from satellite tracking systems. For establishing satellite orbits, both the posi- 
tion and velocity vectors of the satellite must be known. Usually, these two vectors are 
obtained from observations of the satellite made with various types of tracking systems. 
Since no tracking system can be considered completely free of errors, the observations 
obtained from such systems obviously contain errors, of which the most important are: 

(1) Tracking system errors, 

(2) Station position errors, 

(3) Trajectory errors; i.e., satellite orbit errors. 

These three kinds of errors will be discussed in detail here. 

A separate treatment of errors due to atmospheric and ionospheric refraction (Ref- 
erences 1-6) will not be necessary. These errors are assumed to be included in the system 
rms errors, that is, in the errors of the measured quantities such as range, range rate, 
angle, and angular rate. 

To avoid too voluminous a paper, the work has been divided into two separate parts. 
Part I mainly treats the theory; Part II (in preparation) will treat the application of Part I 
by means of high speed computers, and will present concrete examples. 
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Part I is devoted entirely to an error analysis of tracking systems. The systems con- 
sidered are complexes formed by the combination of radars; angle and angular rate systems 
(e.g., interferometers); and range and range rate systems. The errors sought in this analy- 
sis are those resulting from the tracking system itself as well as those due to the uncertainty 
in the systems location. The errors for these tracking systems have been assumed to be 
uncorrelated, on the basis of practical experience which showed that the slight existing cor- 
relation does not appreciably alter the result obtained under that assumption. To permit 
the calculation of these errors it is necessary also to assume knowledge of the orbit or 
trajectory. The latter assumption imposes no restriction on an analysis of this type, since 
the errors considered here are independent of orbit or trajectory errors. It has further 
been assumed that the tracking systems and complexes are capable of determining the posi- 
tion and velocity vectors independently of an assumed orbit. In brief, for some of the dis- 
cussed systems, simultaneity is the minimal constraining condition for the solution. For 
instance at least three range and range rate systems, or one radar system, is required. 

Nevertheless, to facilitate computation, an orbit has been assumed. (Available orbit 
generators permit immediate evaluation of the partial derivatives involved in the computa- 
tions.) Because of the independence of the orbit — at least for a first order approximation - 
approximate orbits or trajectories derived from nominal injection parameters (provided in 
all satellite system operations plans) or from available observational tracking data will be 
adequate for the error analysis. 

The aim of this paper is to obtain equations which permit the calculation of all errors 
resulting from any combination of tracking systems (such as radars, interferometers, or 
range and range rate systems). The coordinate systems used were so chosen that these 
errors could easily be transformed into orbital element. 

Part II will apply the equations derived here assuming simultaneous and non- 
simultaneous observations. Numerical examples will be given illustrating the propagation 
of errors for discrete points and intervals along an orbit. The influence of the number as 
well as the distribution of tracking stations around the globe on the position and velocity 
errors will also be discussed. Treated will also be the minimum number of observations 
necessary for establishing an orbit of a stipulated accuracy with given tracking systems 
and their global distribution. 


A. DEFINITION OF SYMBOLS 

The following symbols will be employed throughout the development: 

<£ Geodetic latitude of the tracking station 
\ Geodetic longitude of the tracking station 
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h Height above geo id of the tracking station 

Equatorial radius of spheroid used to represent the earth, (for Hayford Spheroid 
a = 6378388 meters) 

tf> 

t , 2 Square of eccentricity for spheroid, (for Hayford Spheroid e/ = 0.0067226700223) 
Greenwich Sidereal Time at 0 h Universal Time (U.T.), obtained from the various 
almanacs 

0 Greenwich Sidereal Time at U.T. of observation 
t U.T. of observation 

a Azimuth of the object being tracked 

6 Elevation of the object being tracked 
r Slant range of the object being tracked 
a Semimajor axis of satellite orbit 

;v' Right ascension of satellite 

h Declination of satellite 

C Eccentricity of satellite orbit 

s 

i Inclination of satellite orbit 

s 

n Longitude of ascending node 

^ Argument of perigee 

s 

E Eccentric anomaly 

s 

p Period of revolution of satellite 

s 

N Radius of curvature along Prime Vertical 
p Magnitude of station position vector 

p Magnitude of the radius vector to the satellite in the inertial coordinate system 

B. MATRIX NOTATIONS 

Since a large portion of this paper deals with transformation of coordinate systems, 
it will be useful to present a brief introduction to the transformation theory relevent to the 
material presented herein. A translation indicating a shift in coordinates is represented 
in matrix form as 


X 


Y-S 


(B- 1) 
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or, expanded in components, 


( x i- x 2 ’ x 3 ) = (yi - s i ’ y 2 " v " s 3 ) ■ (B-2) 

A rotation on the other hand indicates an angular shift of two coordinate axes about a 
third axis — the chosen axis of rotation. In matrix notation a rotation is represented as 


Z = Ri (r)Y, (B . 3) 

where y is the argument of the rotation and the index i refers to the axis of rotation in the 
y- coordinate system. For instance, 


I yr) 


1 

0 

0 


0 

cos y 
- s in y 


0 

s in y 
cos y 


(3x3) 


(B-4) 


The rotation illustrated in Equation (B-3) represents a rotation of a Y -coordinate 
system around its yj -axis through an angle y. This similarly holds for rotations around 
the y 2 or y 3 -axes. 


In this paper extensive use is made of rotation matrices; therefore, it is now appro- 
priate to mention the following useful properties of such matrices: 


( 1 ) 

RT ( 7 ) = R;- 

(y) 

( 2 ) 

R ; T ( 7 ) = Rj 

(-y) 

( 3 ) 

Ri (7) RJ (7) 

= I 


(B-5) 


where i = 1, 2, or 3, and I is the identity matrix, which can easily be seen from Equa- 
tion (B-4). 


C. DEFINITIONS OF COORDINATE SYSTEMS 
1. Coordinate Systems Used 

The following coordinate systems are used: 

X Inertial Cartesian coordinate system (Figure 1); 
X S [Xj, x 2 , x 3 ], 

X j -axis directed towards vernal equinox, 
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x 3P y 3 



(GREENWICH) 

Figure 1— Geometric representation of the coordinate systems referred to the earth 

ellipsoid. 

x 2 -axis normal to both x t and x 3 -axes, 

x 3 - axis directed along earth’s axis of rotation. 

Y Earth-centered and earth-fixed coordinate system (Figure 1); 

yj -axis directed towards Greenwich, 

y 2 -axis normal to (y 3 , y 3 ) axes, 

v -axis directed along earth's axis of rotation. 

¥ Coordinate system parallel to the Y-system, centered at the observer. 

Z Local Cartesian coordinate system centered at the observer (Figures 1 and 2); 

z x -axis directed towards local East, 

z -axis directed towards local North, 

2 

2 —axis directed along normal to local horizon plane, 
s Position vector of observer with respect to Y-system (Figure 3 and Appendix A). 



Figure 2 Cross-section of the elliptical assumed earth (Hayford ellipsoid) 

2. Relationships Between Coordinate Systems 

Some relationships between the various coordinate systems discussed should also be 
given: The transformation from the x-system to the Y-system by simple rotation is 

Y " R 3( t9 G ) X ' (C-l) 

where 

& r = 0 r + i?t , 

G G 0 

= 0.26251595 radians /hour 

and 


0 < < 2rr . 


SATELLITE 

S a A" 



Figure 3— Artificial satellite orbit as referred to the local horizon coordination system 


The transformation from the Y -system to the Y -system by translation is 


Y = Y-S , 


(C- 


where 


and 


s ^ S 


(N + h ) cos A cos K 
( N + h ) cos <j> s in K 
|n (l - Cj 2 ) + h ] sin <f 


N 


[l - e * s in 2 4]'' 2 


(C- 


The transformation (double rotation) from the Y-system to the Z-system (Figure 1) 
can be written 
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To simplify the notation, let 



then the transformation of the \-system into the Z -system can be written: 

Z = H R 3 (t? Q ) X - R S , 

or in partitioned matrix form 


Z 


(3xl) 




(C-5) 


(C-6) 


(C-7) 


Next the coordinates of the satellite in its orbital plane U will be transformed into the 
Cartesian coordinates of the inertial system X. (Figure 4) In this paper it is assumed that 
the earth's center of mass and the center of the spheroid representing the earth are coinci- 
dent. This, in turn, implies that the origins of X and V coincide; hence only rotations are 



Figure 4— Orbit in the inertial frame of reference 
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necessary to perform the required transformation. The equation for this (via the Eulerian 
angles) is expressed as follows: 


X = R 3 (-0) Hj ( i ) R 3 (-*) li (C-8) 

where — 

a s (cosE s - ej 

u = a s _e s 2 )‘ 2 sinE s 

and the elements of U actually represent the Cartesian coordinates of a Keplerian ellipse in 
the orbital plane. 

D. MATHEMATICAL MODELS OF TRACKING SYSTEMS AND THEIR ERRORS 

The primary purpose of any tracking system is to determine the position vector r It ) 
and the velocity vector v = dr dt of an object moving in space. The position vector can 
be expressed as 


where r - |r | is the magnitude of r and r° is the unit position vector. Differentiating 
Equation (D-l) directly gives the velocity vector 7 of the object: 


dr 

Tit 



dr° 
+ r dt 


(D-2) 


The four basic tracking systems which exist and which satisfy Equations (D-l) and (D-2) 
will now be described. 


The Radar System 

Herein it is assumed that a radar measures the range r (which in actuality is a com- 
puted value based upon the travel time of the radar wave), the azimuth a, and the elevation 
angle e. From these, the position vector r is fully determined. The time derivatives r, 
i and k are excluded in this treatment of a radar system because of the relatively large 
errors involved. The latter omission, however, should not give the impression that a radar 
system measuring r, a, e , i, a and i cannot be treated by using the equations developed 
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in this paper. For such a case, the equations for one range system and one angle and an- 
gular rate* system should be used. 


The Angle Measuring System 

An angle measuring system is best illustrated by the radio interferometer (Refer- 
ences 7 and 8) which measures the difference in the arrival times of the wave front from 
a distant point source at a pair of receiving antennas separated by a known distance or 
baseline . This radio path difference is measured by comparing the phase angles of the 
signals received at the two antennas. Two such baselines are employed (four antennas)' 
and from the resultant differences in phase angle A <p along both baselines, the direction 
cosines for two components of the unit position vector r° are determined from the following 
relations: 


n - (l - / * _ *2) i -2 


where d ^ is the baseline length (in electrical degrees) directed, for example, in the east- 
west direction, and d 2 is the baseline length (in electrical degrees) directed in the north- 
south direction. From the calculated values (/, m , and n) the azimuth a and the elevation 
angle e, can be derived, or, in the inertial coordinate system, the right ascension and 
the declination 8. 

The unit vector r° is obtained at discrete times t., where i = 1, 2, ••• n; therefore 


dr° AF °( t 1 ) 
dt * At. 


can also be obtained. This means that each observation consists of (a, a, k) at t . For 
a good determination of r and v which satisfies Equations (D-l) and (D-2), simultaneous 
observations from two stations are required to determine the three components of r and the 
three components of v . In actuality this is an over-determined set of observation equations, 


♦Angular rate information, not measured by a radar system, can be calculated numerically (also see the section on smoothing) 
by fitting a polynominal to the measured angular information. 
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since the number of measurements is greater than the number of unknown parameters; that 
is, the least squares condition is satisfied and could be applied for determining r and v. 


The Range Only System 

It is presumed that this system measures only the slant range r in the manner de- 
scribed next for the radar system. Because only one measurement r is made at time 
t. ( i = 1, 2, ••• n) a range system is inadequate for determining the three components of 
the vector 7. To make it possible to determine 7 at t., simultaneous observations must 
be made from at least three tracking stations whose positions are known; it is also pos- 
sible to determine roughly the vector (Reference 6): 

_ . *(0 

v " At. ‘ 

1 


from discrete measurements r(t.). 


The Range Rate System 

A range-rate system measures d/dt |7|. To determine the velocity vector v, at 
least three simultaneous observations are required from three stations in order to deter- 
mine the components of v and satisfy the vector Equation (D-2) (see Reference 15). As in 
the foregoing determination, the quantity "determined", i in this case, is actually a calcu- 
lated quantity; it is obtained from observing the Doppler shift Av Q of a frequency v 0 . For 
a first approximation r = c(Av 0 /v 0 ) where c is the velocity of light (299.7929 x 10 m/sec). 

Advanced tracking systems are combinations of the aforementioned basic systems 
for example, Azusa and Cyclops. In order to treat these various types of tracking systems, 
their mathematical models, derived from the 7 and 7 vectors, are now presented and the 
corresponding variational equations lead finally to the determination of tracking system 
errors. 


1. The Radar System 

A radar system fully determines the satellite position vector 7 whose components in 
the local Cartesian coordinate system are (z, , z 2 , z 3 ). We shall now derive the varia- 
tional equation relating system uncertainties in range, azimuth, and elevation to the com- 
ponents of 7. 
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Since it is convenient, in this case, to 
work in spherical coordinate systems, the 
components of the position vector r in the lo- 
cal coordinate system (shown in Figure 5) can 
be written in matrix notation as 


Z 


(3x1) 


Vz 1 


r sin a cos € 

Z 2 


r cos a cos € 



r sin e 

_ Z 3_ 

(3xl) 

_ _ 


(D-3) 


The variations (8r) of r, which will be 
used later for the error calculations, are de- 
rived from the first order terms of the Taylor 
series as applied to Equation (D-3) and are 
expressed in matrix form by 





Figure 5 — Satellite position with reference to 
local coordinate system 




dz j 

a _± i 

dz i 


r— — 


__ 

da 

de 


5z i 


111 


dz 2 


P'1 

8z 2 . 

- 

dr 

da 

de 


ba 

J Z 3 

(3X1) 

SZ 3 

dz 3 

dz i 


be 



dr 

da 

de 







(3X3) 



By evaluating the coefficient matrix in Equation (D-4) from Equation (D-3) we obtain 



cos e s in a 

r cos e cos a 

- r s in £ s in 

a 


8 r 

3X 1 ) 

cos e cos a 

- r cos e sin a 

- r s in £ cos 

a 


8a 

sin € 

0 

r cos £ 

— 

(3x3) 

be 


or 



cos £ s in a 

- cos a 

-sine sina 


8 r 

S ^( 3X1) 

cos e cos a 

s in a 

~ sin £ cos a 


- r cos e 8a 


sin £ 

0 

COS £ 


r be 





(3x3) 

— 1 


(D-5) 


(D-6) 


where the coefficient matrix of Equation (D-6) is now an orthogonal matrix J which gen- 
erally simplifies the necessary matrix operations. In matrix form Equation (D-6) now 
becomes 


6Z = J D 6K , 


(D-7) 
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where J is the orthogonal matrix in Equation (D-6) which contains only the trigonometric 
functions of 6 and a, and the matrices [) and sk are defined as 


D 


6 K 


1 

0 

0 


0 0 
“ r cos e 0 


0 


r 


—* (3X3) 


Sr 

S a 
be 


The matrix D is introduced only to separate computational values from values based upon 
observational errors such as Sr, be, and ha. The known values of Equation (D-7) are j, D 
(from an assumed orbit), and 6K (from a tracking system). Because of the orthogonality 
property of the J- matrix (J T J I is the identity matrix), it follows from Equation (D-7) 
that 


D (3x3) 6K (3x1) “ J( T 3x 3 ) 5Z (3xl) * (D-8) 

The form of Equation (D-8) will prove to be more desirable in combination with other 
types of tracking system equations to be used later. Only the position vector r is fully 
determined by Equation (D-3). The velocity vector v = r , on the other hand, cannot be 
precisely determined for a radar system because of two factors; (1) The relatively poor 
angular rate data ba, be (Reference 7); and (2) The relatively large uncertainties in the 
incremental range measurements which are to be used for the velocity determination 
(References 9, 10, and 11). 

2. The Angle and Angular Rate Measuring System 

Variational equations reflecting uncertainties in azimuth, azimuth rate, elevation, and 
elevation rate will be developed now for the purpose of relating these uncertainties to the 
components of the r and v vectors of the satellite in the local Cartesian coordinate system. 

The unit position vector r° can be written by using Equation (D-3) with r = l: 



s in a cos eH 

1 

Z 1 


l 


cos a cos e 
_ sin € J 

I'- 

ll 

X 

Z 2 

(3X1) 

m 

n_ 


(D-9) 
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Since n - (1 - Z 2 -ra 2 ) 1 ' 2 is a derived quantity (in Minitrack, Azusa, Mistram, etc.), 
Equation (D-9) should be restated as follows: 


(2xl) 


s in a cos e 
cos a cos e 


(2xl) 


A variation of the above equation is 





l(2xl) 


cos a - sin a 
~ s in a - cos a 


(2x2) L 


COS € da 

sin € be 


J (2X1) 


From Equation (D-9) it is seen that 


Y 

_1 

2 i 


r 

'(2x1) 

_ Z 2_ 


Z “ + Z , 


(2x1) 


and, again, the variation equation can be written*: 




2 

"(i-i 2 ) 

- lm 

~ in 

bm 

'(2x1) 

r 

- lm 

1 - m 2 

-m _ 


-"(2x3) 




b z _ 


£>z„ 


J (3xl) 


By combining Equation (D- 11) with Equation (D-12) to give 

Sa = F(dz j , Sz 2 , ^> z 3 ) ’ 

Se = g(5zj, Sz 2 , Sz 3 ) . 

we obtain the following result: 


(D- 10) 


(D- 11) 


(D-12) 


- r cos e 8a 


- cos a 

sin a 

0 

r be 

(2xl) 

- s in e s in a 

“sine cos a 

cos e 


(2x3) 


8Z 2 


(D- 13) 


J(3xi) 


•The form of the variational equation requires knowledge of r which is usually obtained from nominal computed trajectory 
information with accuracy sufficient for Equation (D-12); for example r _1 or (r + Ar) -J with Ar/r S 10~ 2 does not alter 
Equation (D-12). 
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Comparison of Equation (D-13) with Equation (D-8)* suggests that Equation (D-13) can 
also be obtained if we multiply Equation (D-8) with an operator matrix F a (where the index 
a refers to angular data) which must be of the form: 


F 


\ 2*3) 


0 1 o 

0 0 1 


Expressing the variational Equation (D-13) by using the given parameters 6K in con 
junction with F a , we obtain*: 


F, 

( 2X3) 


’(3X3) 


®*^(3X1) 


F *(2X3) J ^ X3 > 8Z(3X1) 


(D-14) 


from Equation (D-8). 


For simplicity, as well as consistency with Equation (D-8) (which will be useful later), 
we let 


and 


D° - F D 

a (2x3) 


0 - r cos c 0 

0 o 


-1( 2X3) 


J„° s f J T = 

(2X3) 


0 sina ” cos a 

sine cos £ cos a cose sin a 


(2x3) 


(D-15) 


Now Equation (D-14) can be written: 


D„° 6 K , 3 xn 

8 ( 2x3 ) < 3X, > 


J a ° T 6Z (3X1) 

a (2x3) K ' 


(D-16) 


which is the form of our basic variational equation (D-8). 

The operator matrices such as F a (here) and F r (in succeeding sections) were intro- 
duced for notational purposes only. On the other hand this complicates the computational 
scheme owing to redundant computer multiplications. This difficulty will be avoided by 
complementing the generalized equation with an equation more suitable for computer use f 

•Since Equation (D-8) fully determines Si, it becomes out fundamental equation. 
tThis was suggested by Mr. R. Sandifer. 
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From an interferometer type tracking system, both angular and angular rate data are 
obtained with great precision. Therefore, an equation describing angular rates is necessary 
for a complete description of the interferometer system. 

Utilizing the fundamental Equation (D-7) we obtain the following by differentiation with 
respect to time: 


^ (3x1 ) J(3X3) ®(3x3) *( 3 xl) i 


(D-17) 


which is a direct analog to (D-7) but contains the time derivative of the parameters r, a, 
and e instead of their variations. From Equation (D-17) and using J T J = I, we can calculate 


D K - T T 7 

(3x3) (3X1 ) - J ( 3 x3) ^(3X1 ) » 


(D-18) 


which is consistant with the fundamental variational Equation (D-8). 
Varying Equation (D-18) and rearranging the form we obtain 


®(3x3) 6K (3xl) 


J T )Z - (8D)kl xl) + J ( * x3) tt (3xl) 


where 


6K 


(3xl) 


“S t 
fea 

uu 


(3*1) 


> (D-19) 


J 


By substituting for K from Equation (D-18), calculating 8J T , sD, and combining terms, we 
have the following equation: 


D SK = I T S7 4 V T T c7 

( 3 x 3 ) ^(3X1) J(3X3) OZj (3xl) T (3x3) J (3x 3 ) <>/j (3xl) 


(D-20) 


which is the form consistent with Equations (D-8) and (D-14). The matrix V = [v i .] in 
(D-20) now has the following elements: 


v i t = 0 , 


cos a - Zj sinaj 
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- 7^j sin e sin a + z 2 sin 


sec e / . . \ 

— p — ( Zj sin a + z 2 cos a j 


tan e 


l “) 


-pfzjCosesina + z 2 cos £ cos a. + z 3 sin ej 


v., - - v„ 


3 1 ’ 


23 32 


V,, , 


which were obtained by expanding Equation (D-20) into its component form. It can be seen 
from above that matrix V, is skew- symmetric; that is v ;j =-v... 

Equation (D-20) can now be modified in a manner similar to Equations (D-8) and (D-16) 
to contain only the angular parameters a, e, and not r. 

The range r was included in Equation (D-19) because the matrices involved are square 
non-singular ones, which simplifies the matrix algebra necessary to obtain Equation (D-20). 
By multiplying both sides of Equation (D-20) by the matrix F ^ and then using Equa- 
tion (D- 15) we have the variational equation, in matrix form, which relates uncertainties 
in azimuth rate and elevation rate to uncertainties in the components of the vectors r and 
v of the satellite; that is: 


D° 6K 

a (2x3) 


(3x1) 


J° 


sZ 


"(2X3) 


(3x1) 


F a (2x3) V (3 x3) J(3x3) 6Z (3x1) 


(D-21) 


where 


D°6K 

s 


- r 5a cos e 
r 5 k 


J ( 2x1 ) 


It should be noted that Equation (D-21) is consistent with Equations (D-8) and (D-16). 
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3. The Range Measuring System 

A variational equation relating the uncertainty in range to the uncertainties of the 
components of the satellite position vector will now be derived. Since |'r| 2 - z T z, the re- 
sulting variational equation is 


6r O*i) r Z ( T 1X3) 6Z (3xi) Z (°!X3) 5Z (3xl) ‘ (D-22) 

The variation &r of r can also be expressed (as for the angle and angular rate system) by 
introducing another operator matrix F r (the subscript r stands for range) which permits to 
use the fundamental Equation (D-8): 


6 r 


0 * 1 ) 


F r 5K , 3 x 

(1x3) 


D, 


(1*3) 


(3X3) 6K (3xl ) 


6K 


(1*3) 


(3xl) 


(D-23) 


where (see page 16) 


F = (1 0 0] 


(1*3) 


Similarly, the matrix Z ( 0 j T x3) in Equation (D-9) becomes, in terms of the J (3x3) 
Equation (D-7), 


r\ 1 T 

7 ° r r I T ^ * 0 . 

0*3) f (lx3) J (3*3) J (l*3) (D-24) 

Equation (D-22) can be restated consistent with Equations (D-8) and (D-16) by using 
Equations (D-23) and (D-24): 


( 1 * 3 ) 


8K 


(3xl) 


J r° 6Z., X ., 
(1X3) ( 3xl ) 


(D-25) 


4. The Range Rate System 

Next, the variational equation relating the uncertainty in range rate to the components 
of the satellite position and velocity vectors will be derived in matrix notation. Instead of 
considering the variation of Equation (D-25) we shall consider its time derivative instead, 
so that Equation (D-24) is combined with Equation (D-25) and replaced by 


Fr (lx 3 ) K < 3x1 > 


(1*3) 


i T 7 

J (3X3) *"(3X1) • 


(D-26) 
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The resulting variational equation can be obtained immediately by application of the opera- 
tor matrix F from Equation (D-23) to Equation (D-20): 

r (l x 3) 


D r ° 

0*3) 


SK 


(3X1) 


Jr° 


6Z 


(1x3) 


(3X1) 


+ F r 


(1*3) 


V 1 T 

T (3x3) J (3x3) 


s Z, 


’(3X1) 


(D-27) 


Equation (D-27) represents the desired form of the range-rate variation. 

Considering only those errors (&r, hi- , 8a, 8a, Se, be) pertinent to the tracking systems 
discussed above, we are able to calculate the error components in the local Cartesian co- 
ordinate system Z for each individual system. 

As an analogous problem let us consider a general set of linear equations of observa- 
tions which can be represented in the following matrix form: 


Y 


(nxl) 


X X 

/1 (nXm ) (m xl ) 


(D-28) 


where n is the number of observations, n> m, and m is the number of degrees of freedom 
in the system. It is well known that Equation (D-28) can only be solved by the Method of 
Least Squares (Appendix B); and to solve it for X we proceed in the following manner: 

A T Y 
X 


A T AX 

[a t a]-‘a t v . 


(D-29) 


5a. Error Equations for Tracking Systems in the 
Local Coordinate System Z 

From Equation (D-28), when applied to Equations (D-5) through (D-27), we are able to 
write the equations pertinent to the various types of tracking systems by assuming the fol- 
lowing station complex: 

(1) For a Radars * the error equation (Equation D-8) becomes: 

n sk - I t 57 (D-30) 

U (3ax3) 0f \3xl) J(3ax3) ° #J (3xl) * 

where a >l insures over-determination in position. 


•Here a represents a row index and not an angular quantity. 
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(2) For b Angle and Angular Rate measuring systems, Equation (D-16) becomes 


D a ° «, J>n = J° Si ; 

(2/3x3) ( 3x1 ) (2/3x 3) ( 3xl > 


(D - 31) 


and Equation (D-21) is 


„ 2 * * * 6 * * * K/3xn " J a ° 6Z ^ F a V J T 6Z 

(2/3x3) < 3 l ) (2/3x3) (2/3x3) (3x1) 


(D-32) 


where p>2. 

(3) For y Range Measuring Systems, Equation (D-25) becomes 

F r 5 K (3 x 1) " j/ 6Z (3xl) > 

(rx3) 13 (yx 3 ) ( 3x1 ) 

where y > 3 . 

(4) For a Range Rate Systems, Equation (D-27) becomes 


F r 6, <r3*n 

r (5x3) (3X1) 


r (Sx 3 ) V < 3x3 ) 3 (3*3 ) ® Z (3XI) + J r° (Sx3) 6Z (3 x 1 ) 


(D-33) 


(D-34) 


where 5 >6 . 

The foregoing variational equations relate the uncertainties in the parameters meas- 
ured by a particular tracking system to the uncertainties in the satellite position and ve- 
locity vectors. There are two existing ways which can be used to evaluate the unknowns 
in these equations — in this case, the uncertainties associated with the satellite r and v 
vector components: 

(1) The determination of the unknowns (Sz., Sz.) by simultaneous solution implies that 
the number of variational equations is equal to the number of unknown parameters; for ex- 
ample, consider the radar Equation (D-30) for a = 1 . This means that an exact determina- 
tion of the uncertainties in the position vector components of the satellite is made without 
the method of least squares being required. 

(2) The determination of the unknowns (Sz., 5z.) with the method of least squares 

(Appendix B) implies that the number of equations exceeds the number of unknowns to 

be determined; that is, the system is over-determined. It is possible to meet the condi- 

tion of over-determination in two ways: (1) By using n different tracking systems making 

one measurement; or (2) by using one tracking system making n different measurements 

where n is greater than the number of unknowns. From a practical point of view, however, 

the latter is not desirable because an orbit derived in that manner would be poorly deter- 

mined in orientation and shape (Reference 12). 
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As an example of an over-determined system, we shall again consider Equation (D-30) 
for a> l This type of solution for the unknown parameters reduces the uncertainties since 
redundancy is a well known condition for enhancing the precision of measurements. Thus 
an over-determined system of variational equations is always more desirable than one in 
which the number of equations equals the number of unknown parameters. 

All the mathematical models under Sections D1 through D4 were presented in the 
local coordinate system Z. This follows logically from the fact that all measurements are 
made locally. 

Since not all the stations are in close proximity, their station location uncertainties 
6S xl are reflected in the satellite position and velocity errors. The value «S (3xl) , on 
the other hand, is referred to the geocentric coordinate system Y, (as shown in Figures 1 
and 2). Thus it would seem logical to use the Y-coordinate system in the entire analysis. 
This is not true, however, since the position and velocity errors depend upon the earth’s 
rotation (i.e., the station location changes with the earth’s rotation). Therefore all equa- 
tions will be written in the inertial coordinate system X. An additional advantage is the 
ease with which these errors can be transformed into those of the orbital elements if de- 
sired (Part II of this analysis will treat this point in more detail). 


5b. Error Equations for Tracking Systems in the Inertial 
Coordinate System X 

The variational form SZ (JX|) and the variational form of the time derivative « (3 xi) 
of Equation (C-6) relating the local coordinate system Z to the inertial system X will now 
be used for the transformation of the error equations in Section D1 through D4. 


and 


a e r 


«(3X,) = R (3X3) R 3( e c) (3 x3) 6X < 3 * l > + dt R (3x3 


R 3 g ) 

u 


6X 


(3X3) 


(3xl) 


(D-36) 


The equations in Sections D1 through D4 are written in the inertial coordinate system 
X by using Equations (D-35) and (D-36) and by rearranging terms. For a radar systems, 

these equations are: 


(1) For a radars, 

D (3ax 3 ) 6K (3X1) + [j TR6S ](3ax.) = L (3ax3) 6X (3xl) - 


(D-37) 
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where a > 1 and 


"(3ax3) 


[j t rr 3 K)] 


(3ax3) 


(2) For j3 angle and angular rate systems (interferometers), 
D° 8K /3xn + [F J t R6S1 

a (2ySK3'l < 3xl ) L a j 


( 2@xl ) = L ’\ 20x3) 6X (3X1) 


and 


D J 


(2/3x3) 

where £ 2 and 


**Sc 3X1) + [ 


F . V J T R6s] ( ^ x1) = M a °_ M (3xl) + L.° SX, 


(2/9x3) 


(2/3x3) < 3xl > 


L° 


(2/3x3) 


(3x3) 


M 


(2/8x3) ’ 


<W G J 

r 

d 

]~l 

+ ^1 

^J T R 

~d0^ R 3 (®c) 



M a ° 

(2/3K3) 


[f.m] 


(2/8x3) 


(3) For y range systems, 

*[ F .J ,RSS ] 


= L 0 6X. 3x1 . 

(yxl) ( 7 x 3 ) ( 3xl ) 


where y ^ 3 and 


L° 


(rx3) 


[f,L] 


(rxj ) 


(4) For S range rate systems, 


D 0 6K + 

r ( 5x3 ) ( 3X1 ) 


where S ^ 6 and 


[F r VJ T R6s] = M 0 6X + L r ° 6X , 

L J (8xl) (8x3) ( 3X1 ) (8x3) ( 3xl > 


L° 

(8x3) 


[ p . l 1 


(8x3) ’ 


M r ° = f F M 1 

’(t*3) L r J 


(8x3) 


(D-38) 


(D-39) 


(D-40) 


(D-41) 
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6a. Fundamental and Smoothing Time Intervals 

Since any physical measurement requires time, here we shall assume that a basic time 
interval At. is required to make a single measurement which can be, for instance, a range 
r , angles 'a., £ . ,ora range rate fj. In actuality, the measured value m. already 
represents an average value over the time interval At { ; that is 

n i 72 ~5F~ m * dt ■ (D-42) 

1 

where m * represents an instantaneous "non-measurable" quantity. The values of nr, as 
seen from Equation (D-42) already represent a "smoothed" measurement. 

With a time interval r >>At, modern tracking systems are able to measure a large 
number (i » 1) of quantities n. during r. To reduce this number to a practical level and 
also reduce the errors, the concept of "smoothing" is introduced. 

Now the problem arises of finding the maximum smoothing time, say t = r ; (Figure 6). 



Figure 6— Geometry of the smoothing process. 
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For a large value of r 2 (such as one-tenth of a period for a near earth satellite) the 
observations of y r y(r, f, a, a, e, e) are of course nonlinear. Since the "real" values 
y (p -*°°) of these quantities are never known, we are forced to solve this equation by the 
least squares method, since the number of measurements i far exceeds the number of 
unknowns n. For most practical cases a development of the values y is the polynomial 

p 

y K = £ "W ' ( D - 43 ) 

y=l 

where the number of measurements within t is k 11 l, 2, * •• , i; or, in matrix form, 

Y = r r 

*(ixl) i*P ) ( p xl ) > 


where 


"(ixp) 



(ixp) 


Solving Equation (D-43) by using the least squares method Equation (b-12), (Appen- 
dix B), we obtain for the coefficients a y in matrix form: 


r 

(p*i) 


(c T C)-‘ 

(p xp) 


c 


T 

(pxi) 


Y 


(ixl) ’ 


(D-44) 


where 


r ( pxi) - M (pxl> - 

The measure of "goodness of fit" for a single observation is obtained from Equa- 
tion (b-14) of Appendix b: 


y - y ± v » 


where 


y = T^y, 


tLL 

k =1 7=0 




and 


(D-45) 


K— l 

rA) 7" ( y x _ y) 2 


V 


(D-46) 
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Equation (D-45) can also be utilized for an estimation of the maximum smoothing time 
T . As can be seen from Equation (D-45), the "best” value y for an assumed order P of 
polynomial is a function only of the number of measurements i made during the time in- 
terval r = iAt (where At is the constant basic time interval necessary for making a single 
measurement y). Therefore an arbitrary increase of the time r cannot be made, since the 
values of y, calculated from a series of the order P , would deviate more and more from the 
nominal values (Figure 6). For example, a second order polynomial (p = 2 ) cannot it a 
third order curve over any interval. Thus a criterion for the maximum tolerable smooth- 
ing time for a given polynomial of degree P is 


where i, < i, and i 2 < i 3 etc. Should the above inequality reverse itself, the order of the 
polynomial must be increased (Figure 6). 


6b. Introduction of Smoothing into the Error Equations 


During the course of this error analysis it was assumed that the uncertainties Sr, 8a, 
5e, Sf , &a, and si for the different tracking systems were known and were represented by 
the systems’ uncertainty matrices SK (3x1) and6k (3xl) Equations (D-7) and (D-19). 

The term "smoothing process" means simply that the aforementioned uncertainties are to 
be improved (reduced) in the statistical sense. This means that each value Sr, Sr ■ • • must 
be divided by the square root of the number of measurements i made during the time in- 
terval r. Expressed in matrix form and applied to 6K and sk, the smoothed error equations 

are 


SK 


* 

(3X1) 


^(3x3)®^( 3 xl ) 


and 


sk 


* 

(3x1) 


^(3 *3) SK (3xl) 


(D-48) 


where 


w 

” ( 3 >3) 


1 

fi 1 (3X3) 


and I is the identity matrix. 
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The time of smoothing is then given by* 

r = iAt . (D-49) 

It should be noted that the value of i is finite (uncertainties do not disappear completely) 
since, as was previously mentioned, r and At are finite. 


7a. Error Equations in Position and Velocity Assuming Random Errors 

By using the error equations presented in Section 5b it is now possible to construct 
two generalized equations providing rms-errors in satellite position and velocity, includ- 
ing smoothing. The error equation for position errors is 

A (P*1) = B (p*3) 6X (3xl) ’ (D-50) 

where 


B (3ax3) 5K (3xl) + [j T R 6 ^]( 3 cl 


xl ) 


<P*1) 



SK <‘»> * [f.J 1 ««],„, 


(Pxl ) 


and 


L 


(3^X3) 


with p> 3 , where 


B 


(px3) 


L 0 


(2/3x3) 


L 


(px3) 


p “ (3a + + y) . 

For every case three uncertainty components of the station are always included. Thus one 
range system (y - l), for example, is sufficient to determine the position error, since an 


*In this equation the quantity i is not to be confused with the inclination i of an orbit. 
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orbit is presumed. This statement is academic, however, since this is not a practical way 
to determine spacecraft position errors. Knowledge of the orbit is only necessary for 
evaluating the coefficient matrices D, J and V of sX,6S and 6K, respectively. For an error 
analysis of this type any nominal orbit can be used for evaluating the coefficient matrices, 
since an error in the assumed orbit does not influence the position and velocity error de- 


rived herein; i.e., the errors do not depend on the first order of the "real” position or 
velocity. 

The solution of Equation (D-50) by the least squares method is 


(3x1) (3X3) ®(3xp) A (P*1) 


(D-51) 


which represents, in the inertial coordinate system, the total variations in each component 
of the satellite position vector. 

As was previously mentioned, the uncorrelated errors in the components of the posi- 
tion vector ? are desired (Figure 5). These errors can be obtained from 


H 


(3x3) 


6X, 


3X1) 


6 \ 


T 

(1x3) ’ 


(D-52) 


which contains the correlated terms; H' is formed by assuming that the correlated terms 
are zero, that is, 


h; 


(3X3) 


0 0 

T) 2 0 

*2 

0 V* 


(D-53) 


J (3x3) 


The total rms error in satellite position is then given by (Reference 15). 


L< 


(D-54) 


The position errors determined above can now be used in determining the velocity 
errors, for which the matrix equation is 


(qxi> 


B (qx3) 5X (3X1) * 


(D-55) 
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where 


(q*i) 


<■■*> * t F - v j t 11 «]«... - 

- v j ti '* s U, 


D 0 6 K * ■+ I F 

L r (Sx3) < 3xl > ir 


(qxl ) 


L 0 


(2/3x3) 


(qx3) 


L 0 


r(8X3) J(qx3) 


and* with q > 3 


q { 2fi + b ) 


Solving Equation (D-55), as before, by the least squares method, we have 


6X 


(3xl) 


B ) (3X3) B ^3xq) A (qxl) * 


Thus the velocity errors are obtained: 


H (3X3) 6X (3xl) SXT (ix3) 


And again assuming no correlation in the matrix: 


h: 


(3x3) 


0 o 

r l? 0 

*2 

0 r) ? 


(3x3) 


Therefore the total rms error in satellite velocity is given by 


z>. 


(D-56) 


(D- 57) 


(D-58) 


(D-59) 


If the position error components are not computed from Equations (D-50) , (D-51), and (D-52), then, of course, q >6 (e.g., 
sir range rate stations or three angular rate stations). 
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The general error equations — Equations (D-50) to (D-59) — can be solved from two 
points of view: 

(1) By letting the indices (a, J3, y, or s) represent number of stations which observe 
the satellite; or 

(2) By letting (a, p, y, or 8) represent the total number of observations made by each 
type of tracking station. 

The former approach will require less computer storage, since the evaluation of the 
partial derivatives - see Equation (D-4) - comprising the coefficient matrices of 6K, 6S, 
and eX is kept at a minimum. For example, if the latter method were chosen for a radars, 
30 sets of partial derivatives would have to be evaluated per second instead of the 3 per 
second for the first approach. It was presupposed in this example that the radar sampling 
rate is 1 measurement every 1/10 second. By introducing the concept of smoothing into 
the first approach a result equivalent to the second is obtained. 

In the case of a tracking system not capable of measuring ? and v independently of 
knowledge of the orbit, then the second method is more desirable. For example, consider 
the evaluation of the position and velocity error components of the satellite using one range 
and range rate system. In this case the error is no longer independent of the orbit, since 
knowledge of the orbit is required, and independence of the orbit can be presumed only in 
the case of tracking systems or tracking system complexes (e.g., three range and range- 
rate stations) which are capable of determining r and v independently of the satellite orbit. 
In some cases, therefore, the precision of the orbit must be incorporated in the position 
and velocity error components of the satellite. 


7b. Bias Errors 

Every measurement no matter how carefully it is made, is subject to errors. Experi- 
ence has shown that repeated measurements of the same quantity do not give the same re- 
sult. The two kinds of errors influencing the results are: (1) constant or bias errors 
whose magnitudes depend on parameters such as boresight errors etc. and (2) accidental 
errors referred to here as "uncertainties" which have no fixed cause and hence follow no 
fixed physical law. For example, a measurement of some quantity can be represented in 
mathematical form as follows: 


If m 0 represents the "most probable value" then (An ± 8m) must represent the errors in 
the measurement. Those errors characterized by An are defined as "bias errors", whereas 
the accidental errors or "uncertainties" are those represented by iSm. 
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One possible way to determine the bias errors is to use a precise orbit (determined 
by other means, for instance optical observations) which is to be compared with the track- 
ing systems measurements. The arithmetic average of the resulting residuals will then 
give a good estimate of the bias errors. A more detailed discussion on the determination 
of bias errors combined with numerical examples will be presented as a separate paper. 


8. Comments on Tracking System Uncertainties 

Proper precautions must be exercised in using the error equations derived herein. 

For combinations of tracking systems simultaneously tracking a satellite or a spacecraft, 
care must be taken in selecting the uncertainties (Sr , St, 8 a, 5a, Se, 81) related to the par- 
ticular systems. Consider the case of one, two or three radars used to determine the posi- 
tion vector of a spacecraft. For one radar, the uncertainties ( »r, 8a, Se) are used. In the 
case of three radar stations the position of the spacecraft is fully determined by the three 
range measurements alone (Reference 6). Therefore only 5r i; 5r 2 and 5r 3 need be con- 
sidered, and the angular uncertainties 8a., 8e. can be assumed to be zero since they are 
net required for this case: 8r « ( r Sa cos e ), or 8 r « ( r be. ) . 

In the case of two radars, a very simple "assumption" can be made — multiply the angu- 
lar errors by the factor 1/2. A weighting of the angular errors of 1, 1/2 and 0 seems rea- 
sonable for these three cases. Similar weighting procedures will have to be considered for 
other tracking system combinations. 

In brief, when the solution is over-determined the uncertainties contributing the most 
to the overall error should be eliminated. 
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Appendix a 


Transformation of Cartesian Coordinates into Spheroidal Coordinates 

From Figure 2 (page 7) we can write the equation for the elliptical cross section of 
the earth as (Reference al): 


yl = (a-i) 

a 2 + a ® i 1 _c i 2 ) 


where the semimajor axis is 


a 


V /J 


(a -2) 


and 


y, = 3 ® COS 4-' 


y 3 


a ® sin< £’ • 


(a-3) 


The unit normal to the spheroid in the y 1 y 3 -plane is 




where 


m 



From Figure 2 it follows that 


COS 0 



(a-4) 


(a-5) 
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And by using Equation (a-3) with Equation (a-5) we see that 

r tan 2 c£' 12 

cos <p \ \ + 7 = 1 

L (i - »,’)• 


or, simplifying, 


tan ' - (l - Cj 2 ) tan <p 


We also find from Figure 2 that 


N cos <p - a cos <p’ ; 


and from Equations (a-1), (a-4), and (a-8), we have 


from which it follows that 


/ 1 “ e ^ s i n 2 0 


y j N cos <p , 


y 3 ■ N (l - e^J sin <p 


» j - y j + h cos (p , 


S 3 ~ y 3 + h s in <p 


and thus 


s j - ( N + h ) cos 4> , 


s 3 - ( 1 ~ e j 2 j + h j s i n <p 


Since in Euclidian space 


S = S 1 + S 3 


(a-12) 
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where 


we have 


5j j^cos A i + s in X j ] 


(a-13) 


s 3 - s 3 k t 


S S 


(3X1) 


tN + h ) cos <p cos X 
{ N + h J cos <t> s in X 
N (l ~ e j 2 j + hj s in <p 


(a- 14) 


(3xl) 
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Appendix b 


The Method of Least Squares 

The discussion which follows is a brief exposition on the method of least squares, 
neglecting the conditions of constraint. 

In order to treat the generalized case of the method of least squares, we assume a set 
of ^-observation equations in the following form: 

y, = x 2 , •••, x k ) , 


f 2 ( x i- x 2' X k) • 


(b-1) 


f (x , 

n\ l 


where the equations are subject to the condition that n > k, which implies that the system 
of equations must be over-determined. 

Since any observation is subject to errors, we shall rewrite Equations (b-1) as 


+ £ i = f i( x i 
y 2 + £ 2 = f 2 ( X l 

y n + = f n( x l 


+ Axj , x 2 + Ax 2 , 
+ ^*1 ’ x 2 + ^2 f 

+ Axj , x 2 + Ax 2 , 


, x k + Ax k ) , 
- • • , X k + Ax k ) , 

•••■ x k + ^k) - 


(b-2) 


where the values e. and Ax j , with i = l, 2, ••• , n and j =1,2, ••• , k, represent the errors 
in the observations y i and the corrections to the independent variables x. . 
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In order to determine the quantities of interest Ax., we assume that the relative errors 
are small, i.e., ^ j /x j « l; then Equation (b-2) can be developed in a Taylor’s series by 
using only the linear terms. The following result is obtained: 


JL, df i 

+ e i - f i ( x i ’ x 2 ' ' ' ' • x k ) + / IhT + (higher order terms) , 


j = 1 


or 


where i = 1,2, •••,«. 

From Equation (b-3), we obtain: 


L ax i 

) =1 


Ax. 


and 


^ di 

- L ^ ^ - 


j = i 


^ f 2 


L'X7**> 


(b-3) 


(b-4) 


i = i 


c ^ f 2 

2>*. 


j =1 


(b-5) 


E * "zL ^7'^- 


J = 1 


where Vi (with i - 1, 2, • ■ ■ , n) represents small residuals resulting from the general 
equation (b-3). 
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Gauss (References bl-b5) made an assumption stating that the "best" values of Ax j 
are those for which 


<P 


n 

w v 2 - minimum . 

i i 

i = 1 



(b-6) 


where w. represents a weighting factor. Equation (b-6) is known as the condition for least 
squares. For the following discussion we shall assume that w 4 - 1, i.e., all measurements 
are to contain errors of equal importance. 

In order to simplify the algebraic notation, the symbolic matrix notation will be used. 
Thus Equation (b-5), restated in matrix form, is 


V = & 

(nxl ) v '(nxl ) 


P( nxk ) k x l ) 


(b-7) 


where 



V 1 


e 


Ax, 




e 


Ax 

v 

T ( nxl) 

V 2 
, : 

ii 

X 

2 

7T 

X 

II 

2 


' 

V 

n | 


€ 

m n_ 

( n Xl ) 

_Ax k _ 


and 



df 1 


dx l ’ 

d* 2 ' 

’ dx * 

**2 


af 2 

3x7 • 

3*2 ' 

dx. 

k 




dx 1 ’ 

3x 2 ■ 

c)x t 


The matrices Vand £ca n be considered as ^-component vectors, whereas Axis a k -com- 
ponent vector. 

The Gaussian assumption, therefore, may be considered in terms of vector algebra 
as the inner vector product which in matrix form simply becomes: 


<P 


V T V - minimum 


(b-8) 
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where 


V T - transposed matrix V of Equation (b-7) . 

By introducing Equation (b-7) into Equation (b-8), we have 

0 = € T £ - <? T P AX - a\ t P t £ + aX t P T P AX ^ minimum . (b-9) 


The minimum of Equation (b-9) can be determined by differentiation with respect 
to AX: 


t90 

<9AX 

0 


0 . 


(b-10) 


which gives 


S T p - ax 0 t p t p = o . 


(b- 11) 


But the minimum condition also requires that 3 2 <p/d aX 2 be positive and from Equation (b-11) 
it is found that 


dH 

d AX 2 


X=\ 


0 


P T P 


If P T P is positive definite, 4> t is a minimum. 

Equation (b-11) is called the "normal equation", and is easily solvable by using the 
standard method of multiplying both sides of Equation (b-11) by the inverse matrix of P T P, 
designated as [ P T P ] “ 


AX a 


V (kxl) 


[ P] (kxk)P (kxn) ^(nxi)’ 


(b-12) 


For ease of notation we let 


0 


( k xk ) 


[p t p] 


-1 

(kxk) 


and write Equation (b-12) as 


AX n 


(kxl ) 


ft P T 

” ( k xk ) ‘(kxn) ^(nxl) 


(b-13) 
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i = 1, 2, • • • , n 

j = 1, 2, • • k 

Equation (b-13) gives the "most probable" values, in the Gaussian sense, of AX at X 0 . 

Next, we shall determine the measure of "goodness of fit", or the root-mean-square 
error -q . Since in most cases the observed data (for measurements such as range, range 
rate, angle, angular rate) will not obey the normal distribution law, the usual symbol o is 
replaced by v . 

By substituting Equation (b-11) into Equation (b-9), and having now evaluated aX 0 , we 
determine v 2 in the usual manner 

* = 1FF VTV = 7rM<? T £-£ T PaX o) . (b-14) 


where 

_ pTf 

P( nxk) " ” ^ X iJ (nxk ; 

and 

^( kxk ) ^^j^(kxk) ^(kxk) 


where n > k, n - total number of observations, k - total number of unknowns, and 77 is the 
measure of "goodness of fit" or the rms-error of a single observation. 

One question still remains: how accurate are the calculated values AX 0 ? To answer 
this, we must consider H: 

AX = AX 0 ± H , (b-15) 


where H is the uncertainty associated with each component of aX 0 . 

To calculate the components of II, we must first evaluate the components of aX 0 indi- 
vidually. This can be accomplished by expanding Equation (b-11) as follows: 




j = i 


/ p i 1 e j ^ X 01 P j 2 l + ^*0 2 ^~' I P jl p j2 + + ^ X 0k p jl p jt 

3 = 1 j = i i - 1 J = 1 

n n n n 

J2 p>i€i = ^oi XI Pj2Pj i + Ax ° 2 H Pj ' 2 + + Ax ° k ZI P]2P, ‘ 


j *1 


j =1 


i =1 


> (b— 16) 


L 

i = 1 


jk € ) ^ x 01 P jk P jl + ^ x 02 p jk p j2 + + Ax 0 


/ I 

j = l 


j = 1 


£ 

j = l 
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Now Equation (b-16) can be solved for values Ax oi by using a simple method of elimination. 
For example, if we select element q,, of Equation (b-13) we can determine Ax 0i as follows: 


11 P * i£ i = Ax oi q n7~! p A + "• 4 qn Ax ok T p,,p ik 

1 = 1 J =1 J =1 


n n 

P J2 £ j _ AX 01 Q 12 P j2 P jl 


+ q , _ A.x rt 
1 1 9 n 


j =1 


*12 ‘ " '0 k Pj 2 P j k 

j =1 


n 

/ , P jk 6 j “ ^*01 q ik Pjk P j 1 

J=1 j=I 


n 

+ “■ + q ik P A 

i =1 


> (b-17) 


Summing Equation (b-17) by columns, gives in the first element of the vector aX 0 
Equation (b-13). Furthermore, by definition of the matrix y given in Equation (b-12), 

Q(kxk) (P T P) ( kxk) = I(kxk) • 

The elements in the first row of this matrix multiplication are: 


71 n 

>ii 22 p i 2 i + q i2 22 p >* p " 


n 

22 p i» p i' 


j = i 


22 Pjip * j + q > 2 
j * i 


* n 

22 p " + + q,k 2] P j 2 p j i 

i=i j=l 


= 0 


n n 

22 P) ip i k + q i2 22 p >* p > 


j =i 


j 


k + ••• + q, k 22 = 0 

j =1 


r (b- 18) 


from which we obtain 


01 ^11 


71 n 

22 p ’' e > + q >* p ^ e i + ■■■ + q ik 22 K ipi 

>*» i=i j=i 




(b-19) 


Similar relationships exist for Ax 02 , •••, Ax . 
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By expanding Equation (b-19), we find an equivalent relationship for AX 01 


Ax oi = £ 


k k k 

1 q i i P l i + e 2 q ii P 2i + * * * + € n /, q i 


(b-20) 


For purposes of simplicity we let 


z: 

i = i 


K 

= L q '' p - 


(b-21) 


so that Equation (b-20) may be rewritten: 

AX 01 = a i £ l + a 2 £ 2 + ••• + a r. e . 


(b-22) 


To use the result obtained in Equation (b-18), we multiply Equation (b-21) by the ele- 
ments of the first column of the matrix P ( nxk j : 


p n a i = q u p i 2 i + q i 2 p n p i 2 + + q ik p ii p it 


P,> a , 


q U P 2 2 l + q i2 P 21 P 22 + + q ik P 21 P 2k 


Pn. a „ 


IllP^l + ^ 1 2 P ril p n2 + ••• + ^ik p nl P n 


(b-23) 


By summing the terms of Equations (b-23), we obtain 


n n 71 n 

Y] Pji a j = q ll 7", P) 2 1 + q !2 J~1 p ii p n + + q lk P il 

j-i j=l i = 1 i = 1 


(b-24) 
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Comparing the latter with Equation (b-18), we find that 


L 

i = i 


22 p ^ a < = ° 

i = x l 

f (b-25) 


2] = ° • 

i=l J 

Returning to Equation (b-21), multiplying each equation with a. and expanding the sums, 
we obtain: 


“l 2 = q il P U a i + q i2 P 12 a i + ••' + nikP,!,®! 
a 2 2 = q il P 2i a 2 + q i2 P 22 a 2 + + q ik P ,l, a 2 




> (b-26) 


a 


2 

n 


q, iP « 




^IlcP* 


J 


Summing Equation (b-26), we obtain 


n 



i =1 


n n n 

q n 22 + q « 2 22 p,2a> + + q,t 22 p > ka j 

i j =i j =i 


and comparing Equations (b-27) and (b-25) we immediately see that 


(b-27) 


n 



j = I 


M 1 


(b-28) 


Similar relationships exist for the other main diagonal elements of the matrix 0. 
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In Equation (b-22) we have, corresponding to each «.(with j = 1, 2, •••, n), an tj. such 
that the total rms error r, associated with Ax 0I is 


^ + a n r >n] 


(b-29) 


But by definition of vin Equation (b-14) we see that t,, 2 - v? - • ' • - V* v ; thus 


''hj - * 




j = 1 


(b-30) 


which, with the use of Equation (b-28), becomes 


\ = ±7 W 


1/2 


(b-31) 


And likewise, we have 


v 


h 2 


± 




1/2 

22 


n 




± 


v q 


1/2 

kk 


Returning to Equation (b-15), we now can say in general that 

Ax. = A* oi ± V q// 2 * (b-32) 
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